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FURTHER  REFINEMENT  AND  APPLICATIONS  OF  THE  MIXED 
COMPUTATIONAL  ALGORITHM  FOR  PLANE  ELASTIC  CONTACT  PROBLEMS 


E.  Yogeswaren  and  J.  N.  Reddy^ 

Department  of  Engineering  Science  and  Mechanics 
Virginia  Polytechnic  Institute  and  State  University 
Blacksburg,  V A  24061 

ABSTRACT 

The  mixed  finite  element  scheme  developed  earlier  has  been  further 
studied  in  order  to  improve  the  model  characteristics  and  operation  as 
well  as  to  explore  further  applications  with  regards  to  mechanical 
joints.  One  of  the  important  modifications  that  has  been  implemented  in 
the  present  study  is  the  utilization  of  a  dynamic  as  well  as  a  static 
coefficient  of  friction  for  the  evaluation  of  contacting  surface 
behavior  instead  of  a  single  friction  coefficient  as  used  in  the 
previous  study.  This  has  improved  the  model  behavior  and  involved  a 
substantial  reorganization  of  program  methodology.  Other  refinements 
include  the  incorporation  of  a  modified  solution  technique  that  allows 
the  solution  of  the  indefinite  stiffness  equation  system  which  is  formed 
in  the  first  iteration  of  the  first  load  step,  and  the  usage  of  a  finer 
mesh.  The  new  solution  technique  also  avoids  the  halt  of  execution  of 
the  program  whenever  small  elements  are  introduced  in  the  leading 
diagonal  due  to  contact  loss.  _ 

The  case  of  the  pin-loaded  aluminum  plate  has  been  studied  again  to 
obtain  better  correlation  with  experimental  results;  and  a  pin-loaded 
orthotropic  plate  behavior  as  predicted  by  an  analytical  solution  is 
compared  with  the  present  model  predictions.  The  hybrid  technique  has 
been  used  to  estimate  the  static  and  dynamic  coefficients  of  friction  of 
the  aluminum  pin/plate  system. 

An  elastic-plastic  analysis  based  failure  model  of  pin-loaded 
laminates  is  illustrated  by  examples  which  indicate  when  bearing, 
shearout  and  tensile  failure  occur  in  the  laminate  mechanical  joints. 
Finally,  some  improvements  that  could  be  carried  out  on  the  present 
model  to  predict  severe  deformation  and  fracture  behavior  are  discussed. 


^Graduate  Research  Assistant 

^Clifton  C.  Garvin  Professor  of  Engineering  Science  and  Mechanics 


1.  INTRODUCTION 


The  updated  Lagrangian  formulation  based  on  a  mixed  variatonal 
statement  and  the  associated  finite  element  model  developed  earlier  [1] 
gave  results  in  good  agreement  with  the  analytical  solutions  for  most 
contact  problems  studied  there.  However,  the  numerical  results  obtained 
for  the  pin-loaded  aluminum  plate  showed  poor  correlation  with  the 
experimental  results  of  Joh  (2).  The  poor  correTation  can  be  attributed 
to  both  theoretical  and  experimental  models  used.  It  was  also 
discovered  that  the  need  for  a  better  solution  technique  existed  since 
the  program  operation  was  disrupted  whenever  certain  contact  node  pairs 
were  about  to  lose  contact  thus  creating  very  small  terms  in  the  leading 
diagonal.  The  work  reported  here  is  largely  based  on  the  subsequent 
research  carried  out  to  improve  the  accuracy  of  the  results  by  modifying 
and  refining  the  computational  procedure  developed  earlier. 

Applications  of  the  refined  model  to  the  pin- loaded  plate  problem  and 
other  contact  problems  such  as  a  mechanical  joint  in  filamentary 
composite  laminates  are  presented. 

Factors  responsible  for  the  discrepancy  between  the  theoretical  and 
experimental  results  can  be  many  and  may  include  the  non-ideal 
conditions  under  which  the  experiment  was  conducted.  An  attempt  Is  made 
in  the  present  study  to  see  if  closer  agreement  between  the  theoretical 
and  experimental  results  can  be  achieved  by  refining  the  computational 
model.  Indeed,  the  study  showed  that  It  Is  possible  If  a  dynamic 
coefficient  of  friction  (u^)  and  a  static  coefficient  of  friction  (us) 
are  used  in  the  analysis  Instead  of  a  single  coefficient  of  friction,  as 
was  done  in  the  earlier  work.  This  approach  Is  more  realistic  since 
this  comes  closer  to  the  continuously  varying  friction  coefficient 
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observed  in  practice.  The  basis  for  the  values  of  ud  and  u$  is  provided 
by  the  hybrid  technique  study  conducted  with  the  experimental  values  of 
displacements  along  the  hole  boundary. 

It  was  reported  in  [11  that  the  first  node  of  the  finite  element 
mesh  was  selected  as  to  be  constrained  in  both  directions  which  enables 
the  use  of  a  conventional  banded  solver.  This  is  due  to  the  fact  that 
the  leading  diagonal  terms  In  [K11]  (see  [ 1 ] )  will  be  zero  in  the  first 
iteration  of  the  first  load  increment,  thus  giving  an  indefinite  system 
of  equations.  Although  this  simple  solution  worked  for  problems  with 
carefully  chosen  mesh  numbering,  it  breaks  down  when  leading  diagonal 
terms  become  small  due  to  loss  of  contact  between  pairs  of  nodes. 

Hence,  the  earlier  solution  scheme  was  discarded  in  favor  of  the 
technique  suggested  by  Mirza  [3J. 

In  addition  to  the  above  modifications,  which  are  mainly  remedial 
in  nature,  some  applications  of  the  refined  model  were  also  conducted 
particularly  with  respect  to  failure  in  mechanical  joints  of  composites. 
The  advantages  and  pitfalls  of  the  technique  are  also  outlined  in  the 
section  on  applications. 

Finally,  an  alternative  kinematic  description  is  suggested  for 
severe  deformations  and  cracks  encountered  in  contact  systems.  This 
model  when  incorporated  In  the  present  scheme  can  lead  to  economical 
analyses  of  large  deformation  problems  with  fracture  by  avoiding  the 
remeshing  techniques  normally  used  in  these  cases.  Ho  attempt  is  made 
to  Implement  the  computational  procedure  in  the  present  study. 
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2.  PRESENT  STUOY 

2.1  Contact  Stress  Evaluation 

Evaluation  of  the  kinematical  relations  and  tractions  at  the 


contacting  surface  is  a  significant  factor  in  the  analysis  of  contact 
problems.  One  of  the  main  advantages  accrued  from  mixed  formulation  is 
that  the  normal  and  tangential  surface  tractions  can  be  calculated 
directly  from  the  nodal  stress  values  obtained  In  the  mixed  model.  A 
practical  way  of  treating  contact  surfaces  is  to  assume  that  the  forces 
are  transferred  only  at  the  nodal  points  as  concentrated  forces 
resulting  from  the  Integrated  effect  of  contact  stresses  up  to  and 
including  the  contactor  node.  In  the  earlier  work  [11  a  segment  of 
triplet  nodes  were  considered  on  the  contactor  body,  with  the  node  in 
consideration  as  the  middle  node,  and  by  integrating  the  tractions  on 
the  contactor  surface  the  nodal  forces  were  obtained.  In  the  present 
study  a  more  direct  and  simple  algorithm  is  adopted  to  reduce 
calculations  Involved  and  by  using  a  finer  mesh  it  is  made  certain  that 
the  accuracy  of  the  model  Is  not  compromised. 

In  the  present  analysis  the  contact  status  of  any  contact  segment 
of  the  target  body  containing  a  node  of  the  contactor  Is  decided  by  the 
normal  and  tangential  nodal  forces  at  the  given  node.  If  the  tangential 
force  at  a  node  is  more  than  the  frictional  capacity  then  the  node  is 


under  sliding  contact.  The  limiting  frictional  resistance  that  can  be 
sustained  by  a  stationary  node  K  is  given  by. 


Static 
Frictional 
.  Capacity 


).  ■  (“ 


Static  Coefficient 
of  friction  u# 


)/  Normal  Force  \ 
x  (  Component  I  (2.! 

\  at  node  K  /  . 


However,  if  there  Is  already  slip  occurring  in  the  last  load  increment, 
the  frictional  capacity  is  given  by 
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(2.2) 


Although  the  above  concept  was  discussed  in  [1|,  it  was  left  to  the 
current  study  to  implement  this  into  the  program.  Thus  as  a  first  step, 
the  direct  path  discussed  in  this  section  was  chosen  to  check  the 
usefulness  of  this  notion.  It  should  be  noted  here  that  Sections  3.3.4 
and  3.3.5  of  [1J  are  made  obsolete  by  tne  present  approach.  But, 
implementation  of  u$  and  ud  in  the  integration  process  outlined  in  these 
sections  will  be  carried  out  if  the  present  method  is  found  to  be 
unsatisfactory  in  later  applications. 

Sticking  contact  occurs  if  the  frictional  capacity  as  determined  by 
Eq.  (2.1)  or  Eq.  (2.2)  exceeds  the  tangential  force  at  node  K.  The 
conditions  of  contact  from  the  last  load  increment  determines  whether 
Eq.  (2.1)  or  Eq.  (2.2)  is  used.  Any  free  target  segment  that  comes  in 
contact  will  have  a  sticking  contact  since  the  tangential  force  is  zero 
along  the  free  surface.  This  is  also  true  when  a  segment  re-establishes 
contact  after  a  separation.  However,  when  the  friction  coefficient  is 
very  low  the  frictional  capacity  may  be  so  low  that  sliding  may  occur  at 
initial  contact  itself. 

Sliding  condition  is  brought  about  when  the  nodal  tangential  force 
exceeds  the  frictional  capacity  of  the  segment  given  by  the  appropriate 
equation  (2.1)  or  (2.2).  The  node  is  constrained  to  move  only  In  the 
tangential  direction  and  the  frictional  resistance  opposes  the  relative 
motion  of  the  bodies.  The  dynamic  frictional  resistance  opposing  the 
motion  changes  continuously  as  a  function  of  the  relative  magnitude  of 
the  tangential  and  normal  forces  acting  at  a  given  node.  As  a  first- 
order  approximation,  the  value  of  the  resistance  at  the  beginning  of  the 


4 


load  increment  is  assumed  to  be  opposite  to  the  direction  of  motion.  It 
is  to  be  noted  that  the  global  nodal  forces  are  the  external  forces 
acting  on  every  element  to  balance  the  forces  due  to  the  stresses. 
Frictional  forces  are  externally  applied  at  the  contactor  nodes  in  the 
direction  of  tangential  force  and  are  given  by, 

(External  Frictional  Force*)*  =  x  U(j  x 


(normal  \ 
nodal  I 
force  / 


if  and  only  if 


(2.3) 


Normal  <  Tangential 
Nodal  Force  "  Nodal  Force  * 


The  frictional  force  is  determined  at  all  equilibrium  configurations  and 
applied  along  the  target  surface. 

When  the  coefficient  of  friction  is  very  small  the  frictional 
capacity  of  all  the  segments  under  contact  is  identically  zero  at  all 
times.  The  contactor  nodes  follow  the  fixed  target  surface  and  since 
the  target  surface  need  not  be  parallel  to  the  global  axes,  a  local 
coordinate  system  is  defined  with  constraint  of  movement  along  one  axis. 

Separation  occurs  when  the  reactive  contact  forces  act  in  the 
negative  direction  of  the  unit  outward  normal  to  the  contact  surface. 
However,  if  the  forces  as  evaluated  at  the  end  of  an  equilibrium 
configuration  become  positive,  then  there  is  no  contact  force  betwen  the 
contacting  bodies,  and  the  segment  containing  the  node  has  separated. 
This  node  is  considered  to  be  free  and  once  again  is  a  potential 
contactor  node  and  checked  for  contact  overlap  in  subsequent  load 
Increments. 
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2.2  Modified  Solution  Technique 


The  finite  element  equations  resulting  from  the  mixed  formulation 
can  be  written  as. 


where 


IK11!  (k121 

[K12]T  [K22] 


(2.5) 


1 1  3i|; .  3 ij>  •  3^  -  3X  .  3^i*  3i|i  j  3^  •  3^i . 

K  7  .  =  f  [  t - J-  +  t  ( — -  — * 1  +  — -  — i)  +  t  — -  — dxdy 

J  1  xx  3X  3x  xyvax  ay  ay  ax  '  av  j«  i 


yy  ay  ay 


(2.6) 

Thus  it  can  be  seen  from  Eq.  (2.6),  for  the  first  iteration  of  the  first 

load  increment,  that  (K11]  =  0  because  r  =  t  =  t  =  0,  thus 

xx  xy  yy 

resulting  in  an  indefinite  system  of  equations.  Although  Eq.  (2.5)  can 
be  solved  for  this  condition,  by  first  solving  for  stress  and  then  for 
displacements,  a  more  direct  way  has  been  given  by  Mirza  [3].  He 
suggested  that  premultiplying  both  sides  of  the  equation  by  the 
transpose  of  the  global  stiffness  matrix,  one  can  obtain  a  positive- 
definite  system  of  equations.  For  example,  if  Eq.  (2.5)  can  be  written 
as 

(KIM  =  [P]  (2.6) 

$ 

where 


'(01 

(k12i1 

Hu}  ri 

w 

22  *  = 
(k22jj 

{  x  {P}  = 

|{s} 

({0} 

(K|  = 


then  a  positive  definite  system  Is  given  by 

[K1T[KJ{«}  »  IK1{P}, 

and  the  solution  of  Eq.  (2.7)  gives  the  results  of  Eq.  (2.6). 


(2.7) 


This  technique  has  been  extended  in  the  present  study  by  reforming 


Eqs.  (2.5)  into  the  type  given  by  Eq.  (2.6)  whenever  leading  diagonal 
terms  are  small  and  then  carry-*  ig  out  the  operations  described  above. 

3.  LITERATURE  REVIEW 
3.1  Overview 

The  earlier  work  [1]  contains  and  discusses  a  number  of  references 
on  the  finite  element  analyses  of  contact  problems  and  the  present 
review  is  intended  to  be  with-‘n  the  scope  of  the  present  study.  The 
major  aspects  of  contact  problems  have  been  discussed  in  detail  In  [1] 
and  it  would  be  appropriate  here  to  concentrate  on  how  the  model  can  be 
improved  by  modifications  based  on  an  Eulerian-Lagrangian  formulation  as 
discussed  in  [4-6],  and  to  examine  some  of  the  work  that  has  been 
carried  out  on  the  failure  of  mechanical  joints  in  composites. 

Mixed  formulations  of  contact  problems  can  be  further  enhanced  by 
incorporating  a  form  of  kinematic  description,  addressed  in  [4]  as  an 
Eulerian-Lagrangian  displacement  model.  The  conventional  Lagrangian 
description  uses  a  fixed  reference  material  configuration  to  formulate 
equilibrium  equations  where  a  certain  reference  quantity  of  material  is 
"followed"  throughout  its  physical  behavior.  Conversely,  the  material 
can  be  allowed  to  "flow"  across  a  fixed  spatial  reference  and  physical 
measurements  can  be  referred  to  the  fixed  frame.  The  total  Lagrangian 
description  uses  Initial  configuration  as  the  reference  whereas  the 
updated  Lagrangian  description  refers  to  the  current  deformed  state. 

But  the  General  Lagrangian  description  and  Eulerian-Lagrangian 
description  are  more  flexible  in  that  the  configurations  used  to  measure 
spatial  variables  and  to  define  the  strain  measure  need  not  be  the 
same.  Thus,  this  allows  material  flow  across  the  meshes  thus  enabling 
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large  deformations  in  contact  and  fracture  to  be  captured  within  an 
initial  mesh  rather  than  employing  mesh  redefining  techniques.  The 
application  of  this  technique  is  illustrated  in  (5]  for  contact  problems 
and  in  [6]  for  dynamic  fracture  problems. 

Calculation  of  the  strength  and  failure  mode  of  a  composite 
laminate  containing  a  pin-loaded  hole  requires  the  knowledge  of  the  load 
distribution  inside  the  surface  of  the  hole.  Frequently,  a  cosine  load 
distribution  is  assumed  sirce  such  a  load  distribution  greatly 
simplifies  the  calculations.  The  stresses  inside  the  laminate 
calculated  by  a  cosine  load  distribution  may  differ  significantly  from 
those  which  actually  arise  in  the  structure.  As  a  result,  those  failure 
criteria  which  require  an  accurate  knowledge  of  the  stress  distributions 
near  the  surface  of  the  hole,  will  predict  failure  inaccurately  when 
used  in  conjunction  with  cosine  load  distribution.  Work  on  composite 
bolted  joints  has  been  extensive  with  early  studies  concentrating  on 
empirical  design  and  gradually  progressing  to  analytical  methods  for 
stress  analysis  (see  (7-141)  and  the  search  for  appropriate  failure 
criteria  (see  [ 15-20] ).  Empirical  design  methods  proved  too  costly  and 
time  consuming  since  the  large  number  of  variables  in  joint  design 
require  huge  data  bases  for  each  material  and  class  of  laminates.  Thus 
analytical  techniques  have  been  attractive,  although  this  requires 
vigilance  on  the  part  of  the  designer  towards  the  pitfalls  normally 
encountered. 


* 
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3.2  Eulerian-Lagranqian  Description  for  Severe  Deformation  and 

Fracture 

This  section  reviews  the  updated  Lagrangian  scheme  used  in  the 
present  study,  and  outlines  the  Eulerian-Lagrangian  formulation  for  an 
easier  and  accurate  modelling  of  contact  and  fracture.  However,  the 
latter  scheme  is  not  implemented  on  the  computer  during  the  present 
research. 

The  modified  functional  with  the  Lagrange  multiplier  for  the 
updated  Lagrangian  formulation  is  given  by  (see  [1]), 


By  imposing  stationarity  on  in  Eq.  (3.3),  we  obtain  the  following  two 
apprximate  equilibrium  equations  which  form  the  basis  of  the  mixed 
finite-element  model: 


lri j6(ln1  j)dv  +  lSi j5ui , jdv  *  5(lp)  =  Jv  Lri j6(lei j)dv 
1  1  1  (3.5) 

f  ui,j5(lS1j)dv  =  J  °1Jki  lSkt6(lSij^dv  =  0  *  *3,6) 

V1  V1 
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Approximating  the  displacements  and  stresses  by 


n 

ui(x1,x2)  =  z  u^j.(x1,x2) 
lSi j(xi*x2)  =  lSi j4,k<xl*x2^ 


(3.7) 

(3.8) 


and  substituting  Eqs.  (3.7)  and  (3.8)  into  Eqs.  (3.5)  and  (3.6),  we 
obtain  the  mixed  finite  element  model. 


■V1] 


(K 


12,1 


[K12]T  [K22] 


M 

{s} 


(f;-1 

0 


{F} 

0 


NLl 


(3.9) 


where  the  details  of  the  element  matrices  are  given  in  [11. 

A  similar  approach  with  the  Eulerian-Lagrangian  model  would  give. 


nL  -V  <Sj  +  lS1j><le1j  *  l"ij>JdVR 


3U  . 


3U. 


-  J 


lslMj- j'TTf 

3Xk  3XR 


(3.10) 


where 


3xk 

Jkj  =  TXT  and  J 

J 


Similar  manipulations  as  for  the  updated  Lagrangian  formulation  lead  to, 

(3.11) 


^  Ki  1 1  ^12^  [Ki3] 

[ k21 I  [K-2l  [K'3l 

IK^l  IK'21  [k^3] 

where  {u}  is  the  Lagrangian  incremental  displacement,  {U}  is  the 
Eulerian  incremental  displacement,  and  {S}  is  the  incremental  stress 
vector.  The  total  incremental  displacement  is  given  by 

{UT}  =  (u}  -  {U>  (3.12) 
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The  details  of  the  Eu lerian-Lagrangian  method  and  its  application  in  the 
present  context  can  be  found  in  1 4 J  for  the  contact  model  and  in  [6]  for 
the  fracture  model,  although  References  4  and  6  consider  only  the 
displacement  formulations. 

3.3  Failure  analysis  of  Mechanical  Joints  in  Composites 

Stress  analysis,  a  failure  criterion  and  a  strength  model  are  the 
basic  components  of  a  joint  failure  analysis  (or  a  strength  analysis). 
The  three  basic  failure  modes  associated  with  bolted  joints  in 
composites  are  the  bearing  failure,  the  shearout  failure  and  the  net 
tensile  failure  (see  (15-37)  and  (18-23)).  It  has  been  found  that  the 
shearout  and  the  net  tensile  failure  can  be  adequately  modelled  by  a 
plane  elastic  stress  analysis,  with  a  point  stress  failure  hypothesis 
and  a  macroscopic  failure  model  such  as  the  maximum  stress  or  the 
maximum  strain  criterion. 

Recently,  there  is  a  trend  towards  studies  incorporating  ply-by-ply 
failure  analysis  of  bolted  joints,  in  order  to  assess  the  damage  to 
individual  plies  (see  Reddy  and  Pandey  [24]).  This  requires  some  form 
of  macroscopic  criterion  such  as  the  Tsai-Hill  or  Hoffman's  criterion  to 
predict  failure. 

An  altogether  different  approach  has  been  adopted  towards  failure 
by  Hyer  et  al.  ([7], (8)),  who  have  used  the  maximum  radial  and 
circumferential  stresses  as  indicators  of  the  capacity  of  a  bolted  joint 
in  order  to  study  the  effects  of  pin  elasticity,  pin  fit-interference 
and  friction  on  the  capacity  of  a  joint.  They  concluded  that  all  three 
factors  are  detrimental  to  the  capacity  of  the  bolted  joint. 
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4.  APPLICATIONS 

4.1  A  Pin-Loaded  Aluminum  Plate 

Three  different  meshes  were  used  to  model  the  aluminum  plate.  All 
meshes  took  advantage  of  the  symmetry  of  the  problem  and  modelled  only 
half  the  plate.  Further  simplification  was  carried  out  on  Mesh  A  (Fig. 
1),  where  only  a  quarter  of  the  pin  was  modelled  as  was  discussed  in  the 
earlier  study  [1].  However,  both  Mesr»  B  and  Mesh  C  did  not  adopt  this 
simplification  and  they  modelled  half  the  pin  as  shown.  The  contact 
nodal  locations  of  Mesh  A  are  reiterated  here  for  completeness.  These 
locations  are  the  following  (degrees):  0.0,  1.0,  2.0,  4.0,  6.0,  8.0, 
10.0,  12.5,  15.0,  20.0,  25.0,  30.0,  35.0,  40.0,  45.0,  54.0,  63.0,  72.0, 
90.0,  99.0,  108.0,  117.0,  126.0,  135.0,  144.0,  153.0,  162.0,  171.0  and 
180.0.  In  Mesh  B  these  locations  were  spaced  at  9°  intervals  and  in 
Mesh  C  these  were  at  2.5°  intervals  for  the  first  90°  and  at  9° 
intervals  for  the  second  90°  as  shown  in  Fig.  2  and  3  respectively.  The 
number  of  elements  and  the  nodes  in  the  meshes  are  given  as  follows: 

Mesh  A:  236  elements,  286  nodes. 

Mesh  B:  332  elements,  391  nodes. 

Mesh  C:  452  elements,  540  nodes. 

All  nodes  along  the  line  of  symmetry  were  constrained  to  move  only  in 
the  lengthwise  direction  and  the  center  of  the  pin  (numbered  1)  in  all 
three  meshes  was  constrained  in  both  directions.  The  load  applied  was 
distributed  along  the  shorter  edge  away  from  the  pin,  in  the  lengthwise 
direction.  A  dynamic  coefficient  of  friction  =  0.25  and  a  static 
coefficient  of  friction  =  0.35  were  used  in  the  analysis  with  these 
numerical  values  having  been  estimated  from  the  hybrid  technique. 
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The  load  was  applied  in  14  steps  following  closely  the  experimental 
loading  values.  In  the  load  increasing  phase  the  values  of  load  were 
20,  23,  520,  1240,  1460,  1670  and  1980  lbs.  respectively  with  the 
decreasing  phase  values  being  1800,  1600,  1210,  1070,  535,  210  and  40 
respectively  giving  14  load  steps.  The  number  of  iterations  required 
for  each  load  increment  were  8,  8,  3,  3,  3,  1,  2,  2,  2,  3,  2,  3,  3,  and 
3  for  the  case  illustrated  in  Fig.  6,  where  Mesh  A  was  used  and  results 
are  plotted  for  load  step  4,  where  the  total  load  is  at  a  value  of  1240 
lbs. 

The  case  shown  in  Fig.  5  was  analyzed  with  the  modified 
computational  procedure  discussed  in  Section  2.  The  results  show  a 
closer  agreement  with  the  experimental  results  of  Joh  (2j,  than  the 
results  shown  in  Fig.  4,  which  were  obtained  using  the  original 
procedure  of  [1].  It  is  interesting  to  note  that  a  constant  value  of 
friction  coefficient  of  0.15  gives  a  closer  agreement  in  Fig.  4  than  the 
constant  friction  coefficient  of  0.30,  these  values  being  chosen 
arbitrarily.  It  is  possible  here  to  be  misled  easily  to  conclude  that  a 
better  choice  of  value  of  friction  coefficient  might  be  0.15  than  0.30 
unless  the  full  picture  is  revealed  by  Fig.  5,  which  indicates  that 
better  results  are  given  by  a  dynamic/static  friction  model.  Indeed  it 
is  possible  to  conclude  here  that  yet  better  modelling  can  be  achieved 
by  a  continuous  friction  coefficient  variation  such  as  that  given  by  a 
power  law  although  the  Implementation  of  this  may  be  more  cumbersome. 

It  can  also  be  seen  that  the  angle  of  separation  Is  more  realistic  In 
Fig.  5  (after  the  modification)  than  in  Fig.  4  (before  the 
modification) . 
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The  load  decreasing  phase  results  shown  in  Fig.  6  compare  favorably 

with  the  experimental  results,  despite  some  66 %  difference  in  the  values 

of  t  from  40°  to  75°.  This  can  be  considered  fair  compared  to  the 
r  9 

results  given  by  the  method  before  modification,  which  gave  identical 
results  to  the  load  increasing  phase  without  any  negative  values 
for  t 

re 

4.2  A  Pin-Loaded  Orthotropic  -late 

An  analytical  solution  to  the  prob'em  of  pin-loaded  composite 
laminate  has  been  obtained  by  Hyer  et  al.  ( [ 7 ) , [ 8 | )  based  on  a  complex 
Fourier  series  method  and  a  collocation  technique  which  enforced 
boundary  conditions  at  discrete  locations  around  the  hole  boundary. 
Results  were  obtained  by  this  analytical  method  for  infinite  orthotropic 
plates  loaded  by  pins.  These  results  were  chosen  to  be  compared  with 
the  finite  element  results  since  both  methods  capture  the  idealized 
conditions  of  the  model  to  the  same  degree. 

The  mesh  shown  in  Fig.  7  was  selected  to  idealize  the  infinite 
plate  and  the  pin,  with  426  elements  and  509  nodes.  The  nodes  along  the 
line  of  symmetry  were  constrained  to  move  only  in  the  lengthwise 
direction  in  the  same  way  as  the  nodes  along  the  longer  boundary  whereas 
the  nodes  along  both  shorter  boundaries  were  constrained  to  move  along 
the  y-direction.  Normalized  circumferential ,  radial  and  shear  stresses 
along  the  hole  boundary  were  plotted  against  the  radial  angle  and  found 
to  compare  well  with  the  analytical  results  (Figs.  8  and  9). 
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4.3  Application  of  the  Hybrid  Technique  to  Estimate  Static,  and 
Qynamic  Coefficients  of  Friction 

This  technique  basically  consists  of  applying  the  loading, 
prescribing  the  displacements  from  Moire  analysis  to  the  hole  boundary 
and  prescribing  other  boundary  conditions  as  before.  Thus  only  the 
plate  is  discretized  for  this  analysis,  without  the  pin,  and  the  mesh 
used  is  shown  in  Fig.  10.  The  stresses  along  the  boundary  are  obtained 
from  the  analysis  and  are  shown  in  Figs.  11  and  12  p'otted  against  the 
radial  angle.  Thus  the  shear  stress  xr0  and  the  radial  stress  rrr  show 
remarkable  resemblence  to  the  stresses  given  by  Joh  [21  at  a  load  level 
of  1840  lbs. 

In  order  to  assess  the  friction  coefficient  values,  to  be  used  for 

the  regular  finite  element  aQalysis,  t  /<r  ratio  was  plotted  against 

the  radial  angle  for  load  increasing  and  decreasing  phases  and  the 

result  is  given  in  Fig.  13.  In  the  load  Increasing  phase  nearly  all 

contact  is  associated  with  slip  and  the  maximum  ratio  of  t  /t  cannot 

rQ  rr 

exceed  the  dynamic  friction  coefficient  and  thus  u  .  *  0.25  is  a  rational 

a 

choice.  The  negative  ratio  is  maximum  between  55°  and  60°,  where  the 
transition  from  slip  to  stick  occurs  and  thus  u$  *  0.35  is  chosen  to  be 
the  static  friction  coefficient. 

4.4  Analysis  of  Failure  In  Mechanical  Joints 

A  composite  laminate  (0o/±45l,/90°)s  with  laminae  of  the  following 
properties  has  been  used  for  these  studies: 


Ej  *  19.1  x  103  ks 1 

XT  «  229.4  ksl 

T  -  17.3  ksl 

E2  *  2.0  x  103  ksl 

Yt  •  10.1  ksl 

G12  -  0.9  x  103  ksl 

Xc  «  252.1  ksl 

*  0.3 

Yc  «  32.0  ksl 
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It  has  been  established  that  certain  configurations  favor  certain  modes 
of  failure  (181.  This  fact  has  been  used  in  determining  the  plate 
configurations  for  studying  bearing  failure  (Mesh  F,  Fig.  14),  shearout 
failure  (Mesh  G,  Fig.  15)  and  tensile  failure  (Mesh  H,  Fig.  16). 

The  analysis  has  been  carried  out  and  normalized  radial  and 
circumferential  strain  curves  have  been  produced  for  each  failure 
mode.  Failure  is  indicated  by  the  increased  strains  given  by  a 
nonlinear  model  Incorporating  Tsai -Hi  1 1  criterion  as  compared  to  a 
linear  elastic  model.  Results  are  shown  in  Figure  17  for  bearing 
failure.  Fig.  18  for  shearout  failure  and  Fig.  19  for  tensile  failure. 
The  Tsai -H i 11  criterion  used  in  this  study  is. 


+  (/)2 


,%2  ,1  1  1  N  ,1*1  U 

v  ^  ’  ~2  +  *2  ~  ^)ol°2  "  +  ~2  ~  3 

X  V  2  Y  Z  X 

*  1  In  *  ,°V  *  *  ,°6*2  i 

+  ^  ■  ^2)oi°3 +  (r}  +  (r}  +  (r)  -  1 


where  X,  Y  are  either  compressive  (Xc,  Yc)  or  tensile  (Xy,  Yy)  strengths 
and  T  is  the  shear  strength  in  the  xy  plane. 


5.  SUMMARY  AND  CONCLUSIONS 

The  mixed  finite  element  model  developed  in  |1]  has  been  modified 
by  incorporating  a  realistic  interface  friction  conditions  and  a 
solution  procedure.  A  dynamic  as  well  as  a  static  friction  coefficient 
were  used  to  analyze  a  pin-loaded  plate  problem  for  which  experimental 
results  are  available.  The  new  solution  algorithm  not  only  provides 
flexibility  in  numbering  the  nodes  but  also  avoids  the  halt  of 
computation  due  to  the  appearance  of  small  terms  on  the  leading 
diagonal,  during  the  analysis.  An  accurate  contact  stress  analysis  Is 
essential  in  some  applications  such  as  the  study  of  failure  in  bolted 
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joints  of  laminated  composites  and  some  example  problems  have  been 
studied  in  this  area. 


It  is  possible  to  extend  the  capabilities  of  the  present  model  even 
further  by  adopting  an  Eulerian-Lagrangian  formulation.  This  will 
simplify  the  analysis  of  problems,  such  as  cracks  emanating  from  bolted 
joints  by  several  orders  of  magnitude. 
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Mesh  C  for  the  pin- loaded  aluminum  plate 


Radial  and  Shear  S tresses  (ksi) 


Figure  4.  Comparison  of  the  experimental  and  finite  element 

results  for  load-increasing  phase  at  load  level  1240  lbs. 
(the  FEM  scheme  used  is  that  originally  developed  in  [1]). 
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Figure  5.  Comparison  of  the  experimental  and  finite  element 
results  (the  FEM  scheme  used  is  the  modified 
version  of  the  scheme  developed  in  [1];  Load 
increasing  phase  at  load  level  1240  lbs.). 
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Shear  and  Radial  stresses  (ksi) 


Radial  Angle  (in  degrees) 


Figure  6.  Comparison  of  the  experimental  and  finite  element 
results  (the  FEM  scheme  is  the  modified  version) 
for  the  load  decreasing  phase  at  load  level  1210  lbs. 
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Figure  7.  Mesh  D  for  the  pin-loaded  orthotropic  plate. 


Normalized  stresses 


Figure  8.  Comparison  of  the  analytical  and  finite  element 
results  for  the  pin-loaded  orthotropic  plate. 

The  stress  is  normalized  with  respect  to  the  average 
bearing  stress  (81.9  ksi). 
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Normalized  shear  stress 


Figure  9.  Comparison  of  the  analytical  and  finite  element  solutions 
for  the  pin-loaded  orthotropic  plate  (stress  Is  normalized 
with  respect  to  the  average  bearing  stress,  81.9  ksl.). 
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Comparison  of  the  shear  stress  distributions 
obtained  In  the  experiment  and  hybrid  analysis 
of  the  pin-loaded  aluminum  plate  (load  Increasing 
phase  at  load  level  1840  lbs). 


Figure  12.  Comparison  of  the  radial  stress  distributions 
obtained  in  the  experiment  and  hybrid  analysis 
of  the  pin-loaded  aluminum  plate  (load*1840  lbs). 


Radial  Angle  (in  degrees) 


Figure  13.  Variation  of  the  ratio  of  shear  stress- to-radlal 
stress  for  loed  Increasing  and  load  decreasing 
phases  of  the  pin-loaded  aluminum  plate  (results 
obtained  using  the  hybrid  technique). 
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Figure  16.  Finite  element  mesh  H  used  in  tensile  failure  of  composite  joints. 


Radial  Angle  (in  degrees) 


Figure  17.  Radial  and  cicumferential  strains  in  an  orthotropic 

plate  (the  linear  and  nonlinear  models  show  separation 
indicating  bearing  failure  from  e=0°  to  e=35°) . 


Normalized  strains 


Figure  18.  Radial  and  circumferential  strains  in  an  orthotropic 

plate  (  bearing  failure:  0=0°  to  35°;  shearout  failure 
9=  60°  to  80°) . 
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Normalized  strains 
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Radial  Angle  (in  degrees) 

Figure  19.  Radial  and  Circumferential  strains  in  an  orthotropic 
plate  (  tensile  failure  :  9  =  75°  to  90°). 
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The  mixed  finite  element  scheme  developed  earlier  has  been  further 
studied  In  order  to  improve  the  model  characteristics  and  operation  as 
well  as  to  explore  further  applications  with  regards  to  mechanical 
joints.  One  of  the  Important  modifications  that  has  been  implemented  In 
the  present  study  Is  the  utilization  of  a  dynamic  as  well  as  a  static 
coefficient  of  friction  for  the  evaluation  of  contacting  surface 
behavior  Instead  of  a  single  friction  coefficient  as  used  in  the 
previous  study.  This  has  Improved  the  model  behavior  and  Involved  a 
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substantial  reorganization  of  program  methodology.  Other  refinements 
include  the  incorporation  of  a  modified  solution  technique  that  allows 
the  solution  of  the  indefinite  stiffness  equation  system  which  is  formed 
in  the  first  iteration  of  the  first  load  step,  and  the  usage  of  a  finer 
mesh.  The  new  solution  technique  also  avoids  the  halt  of  execution  of 
the  program  whenever  small  elements  are  introduced  in  the  leading 
diagonal  due  to  contact  loss. 

The  case  of  the  pin-loaded  aluminum  plate  has  been  studied  again  to 
obtain  better  correlation  with  experimental  results;  and  a  pin-loaded 
orthotropic  plate  behavior  as  predicted  by  an  analytical  solution  is 
compared  with  the  present  model  predictions.  The  hybrid  technique  has 
been  used  to  estimate  the  static  and  dynamic  coefficients  of  friction  of 
the  aluminum  pin/plate  system. 

An  elastic-plastic  analysis  based  failure  model  of  pin-loaded 
laminates  is  illustrated  by  examples  which  indicate  when  bearing, 
shearout  and  tensile  failure  occur  in  the  laminate  mechanical  joints. 
Finally,  some  improvements  that  could  be  carried  out  on  the  present 
model  to  predict  severe  deformation  and  fracture  behavior  are  discussed. 
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